Parametric Excitation of a ID Gas in Integrable and Nonintegrable Cases 
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We study the response of a highly excited ID gas with pointlike interactions to a periodic modu- 
lation of the coupling constant. We calculate the corresponding dynamic structure factors and show 
that their low-frequency behavior differs dramatically for integrable and nonintegrable models. Non- 
integrable systems are sensitive to excitations with frequencies as low as the mean level spacing, 
whereas much higher frequencies are required to excite an integrable system. This effect can be 
used as a probe of integrability for mesoscopic ID systems and can be observed experimentally by 
measuring the heating rate of a parametrically excited gas. 
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PACS numbers: 67.85.-d,05.45.-a 

The field of ultracold gases has progressed enormously 
toward obtaining quantum systems with desired densities 
and types of constituent atoms, trapping geometries, and 
well controlled interparticle interactions In particu- 
lar, by using either an optical potential or large magnetic 
field gradients, one can confine the motion of atoms to 
one dimension and create interacting ID gases of bosons 
and fermions 0, H|, which can be described by the 
integrable Lieb-Liniger [j| and Yang-Gaudin 0, @] mod- 
els respectively. The purity and isolation of such sys- 
tems from the environment makes them ideal candidates 
for studies of fundamental differences between integrable 
and nonintegrable many-body dynamics. A pioneering 
experiment on this subject has been performed recently 
by Kinoshita and co-workers @. They have shown that 
a ID Bose gas initially prepared in a highly excited state 
does not equilibrate in the lifetime of the experiment, 
whereas essentially the same system with a weaker ID 
confinement thcrmalizcs much faster. 

How do we decide whether a system is integrable or 
not? Let us put aside strict mathematical definitions of 
quantum integrability Q and look at the problem phe- 
nomenologically. What measurement should we perform 
on a system in order to conclude on its integrability? 
The field of quantum chaos suggests to look at its spec- 
tral statistics [llj ■ If energy levels are not correlated [the 
nearest neighbor spacing distribution (NNSD) is Poisso- 
nian] , we are dealing with an integrable or regular system 



llj . If, in contrast, levels repel each other, the system is 



not integrable [12|, ll3[ . Spectral properties were exten- 
sively studied for various systems I3-13|, and, in par- 
ticular, for strongly correlated condensed matter models 
(see [14j for early work) . 

Another signature of integrability is the localization 
of eigenstates of a re gula r system in a certain physi- 
cally meaningful basis lfj, which suggests the dynami- 
cal probe of integrability: if one creates an excited initial 
state localized in this basis, it will stay localized during 
the temporal evolution. In fact, the absence of thermal- 



ization in the experiment [8[ can be regarded as a conse- 
quence of the localization of the Lieb-Liniger eigenstates 
in momentum space. 

In this Letter we compare responses of highly excited 
integrable and nonintegrable systems to an external time- 
dependent perturbation. We explore the idea that a per- 
turbation localized in the same space as the eigenstates 
of the integrable system probes its local density of states, 
whereas in the nonintegrable case the states are delocal- 
ized, and the perturbation, no matter localized or not, 
couples all of them. Considering two ID models on a 
ring we demonstrate that integrable systems can be much 
more stable with respect to slow variation of their Hamil- 
tonian than nonintegrable ones. Namely, we consider the 
model of a single mobile impurity in a Fermi gas and the 
Lieb-Liniger model, and study their response to a peri- 
odic modulation of the coupling constant. This pertur- 
bation is localized in the many-body momentum space as 
it only changes the relative momentum of an atom pair. 
We show that the nonintegrable system is sensitive to ex- 
citations with frequencies as low as the many-body mean 
level spacing, which is exponentially small, whereas the 
threshold frequency in the integrable case is much larger 
and scales polynomially with the system size. 

Consider N atoms with short range interactions in a 
quasi-lD ring-shaped trap of circumference L = 1. If the 
atomic kinetic energies are smaller than the level spacing 
in the direction of tight confinement, the system can be 
envisioned as a ID gas on a ring with the Hamiltonian 
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+ Ya ij 6(x i -x j ), 



i<j 



(1) 



where rrij and the masses and coordinates of the 

atoms. The ID coupling constants g % i depend on the pa- 
rameters of the 3D interatomic interactions as well as on 
the strength of the tight confinement [l6|. Accordingly, 
by changing these parameters in time one can study the 
response of the effective ID system to variations of g 1 ^ . 
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Assume that the system is initially prepared in an 
cigenstate of the Hamiltonian ([1]) with eigenvalue e v and 
eigenfunction ip„, and consider the weak periodic mod- 
ulation g %3 \t) = g lJ + 28 g 13 ' cos{u>t). Then, in the linear 
response regime the probability to remain in the state 
v decreases with the rate fl = 2ir [S(e v ,u) + S(e u , —oj)}, 
where the dynamic structure factor equals 



e v )\(^m v }\ 2 , (2) 



and where F = 2j<j &g lJ &( x i ~ x j)- Exciting an ensem- 
ble of systems (JTJ) leads to a diffusion of the population in 
energy space with diffusion constant 2ttuj 2 S(e 1 uj) result- 
ing in detectable changes of the total energy and entropy. 

The asymptotic behavior of S(e, uj) at small uj gives the 
dissipative part of the response of the system to a slow 
variation of its Hamiltonian and, therefore, measures the 
degree at which this variation can be assumed adiabatic. 
For complex systems it is believed that statistical prop- 
erties of eigenstates, eigenvalues, and matrix elements of 
a perturbation are well described by the random matrix 
theory [ 1 7| . If both H and F were independent random 
matrices drawn from the Gaussian Orthogonal Ensem- 
ble, the average of S{e v ,ui) over an energy interval larger 
than the mean level spacing D(e) would be independent 
of e and uj 18|. Here we show that this low- frequency 
behavior strongly depends on whether we consider inte- 
grablc or nonintcgrable systems. 

Let us consider the model of a single mobile impurity 
interacting with a gas of identical ideal fermions. In this 
case the parameters of (fTJ) are m\ = ... = tun-i = 1, 
?n jv = AI, and g lN = g. It is convenient to work in 
momentum space introducing the Fourier transform 

P1,...,PN 

where all the momenta are integer multiples of 27r. The 
total momentum Q is conserved and keeping in mind 
that pn = Q — Sili 1 Pi we omit the argument p^ in the 
wavefunction ip, which now becomes antisymmetric in all 
of its arguments. The Schrodingcr equation then reads 
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iV-l 



2m 



(3) 



Let us introduce an auxiliary function 

a(pi,...,p N - 2 ) = -V] ip(pi,-,PN-i) (4) 

and using the antisymmetry of ip rewrite the rhs of 
Eq. (0) in the form g [a — Si=i &{■■-, Pi-i,PN-i,Pi+i, ■■■)] 
(Hereafter, for normally ordered arguments we use the 
shortcut a = a(p%, ...,pjv-2))- We then solve Eq. ^ 



with respect to ip and substitute the result into the defi- 
nition of a ((4]) obtaining the equation 
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2p,g 



(1/2/t) sinK 
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N-2 , 
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,Pi-i,Q,Pi+i, ■ 



Li=i Pi) 



(5) 
and 



where fx = M/(M + 1), f = {fJ,/M)(Q 

K 2 = 2^E-^'L' 1 2 P 2 l -Me- 

By solving Eq. ([5]) with respect to E and a we de- 
termine the eigenenergies e v and the "reduced" eigen- 
functions a v . We then calculate the dynamic structure 
factor @ for F = Sg^ 1 ^ 1 8{xi — xpj) by using the 
relation (ip v \6(xi - x n )\i/j v ) = J2 Pl ,..., PN _ 2 a v a >'- Notc 
that Eqs. §5§ and ([5]) conserve parity (simultaneous sign 
change of all pi). The corresponding even and odd ex- 
citation branches are not coupled by F. They have the 
same density of states and contribute equally to Eq. (f2]). 

We have performed an extensive numerical analysis of 
this model for N = 3, 4 in a wide range of energies, cou- 
pling constants, and for different M. We accurately cal- 
culate up to I0 4 excited levels. In the integrable case, 
71/ = 1, we determine e v and ip v from the known Bcthe- 
ansatz solution [lj| and check that both approaches give 
the same result. In Fig. [T] we plot S(s,us) for four cases: 
The upper panels stand for the nonintcgrable case with 
N = 3 (left) and N = 4 (right) with M = m 87Rb /m4 K- 
The lower panels show the integrable case. 

Dotted lines filled vertically to the a;-axis present the 
average of S(e u ,uj)^/e/5g 2 over several hundreds of states 
with odd parity and zero total momentum in the inter- 
val 0.95e < e v < 1.05e. All four presented cases cor- 
respond to the same value of the interaction parameter 

7 = (k 2 re i) a i ~ °- 5 ' whcrc (O = - I) is the 

average square of the relative two-body momentum, and 
o.i = 1/fJ.g is the ID scattering length. We check that 
simultaneous variation of g and e does not change the 
quantity S(e, ui/y/e)-^e as long as 7 stays constant. The 
lower right panel shows three dotted lines obtained for 
e =6.9, 9.1, and 14.2 [xlO 3 ] with the same 7. With 
rescaled cc-axes they collapse to a single curve. The noise 
is due to the averaging over finite number of states and is 
uncorrelated between different curves. We measure u in 
units of the mean level spacing D(e), so the labelling of 
the horizontal axis holds only for the curve e = 9.1 x 10 3 . 

In all our calculations (not only in Fig. [T]) we observe 
that the low-frequency behavior is universal: In nonin- 
tcgrable cases we always see that S(uj) tends to a finite 
constant when uj ~ D, which is consistent with the ran- 
dom matrix theory and with the fact that all states are 
coupled by the perturbation. In contrast, all integrable 
cases are always characterized by a strong suppression of 
S(uj) already for quite large uj/D, which means that the 
perturbation does not couple states with close energies. 
This is the main result of our numerical experiment. 

In order to understand this behavior and the peaks at 
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Figure 1: (color online) Dotted lines present dynamic structure factors for nonintegrable (upper panel) and integrable cases 
(lower panel). Red solid and green dashed lines are the result of the binary approximation (see text). Insets show the nearest 
neighbor spacing distributions for the ensembles of energy levels used to calculate S in the corresponding parent graphs. Red 
dashed and blue dash-dotted lines stand, respectively, for the Poissonian and Wigner-Dyson distributions. 



higher frequencies we use the cluster expansion up to bi- 
nary terms, i.e. we assume that only pairs of atoms can 
be excited at a time, the remaining particles being non- 
interacting spectators. The red solid lines in Fig. [T] cor- 
respond to the quantity S(e,cd) = (N — l)D(e) f Q pie — 
E)S2(E,u})p2(E)dE, where the prefactor N — 1 is the 
number of interacting pairs, S2(E,ui) and P2{E) are, re- 
spectively, the dynamic structure factor and the density 
of states for an atom pair, and p(e — E) is the ideal-gas 
density of states for the remaining atoms. 

For highly excited two-body states with energies w E 
and center-of-mass momenta Q Eq. ([5]) gives two exci- 
tation branches: K nj ± = 27r(n — 1/4) + <f> ± A, where 
6 = — arctan(aifi), n 2 = 2pE — A/£ 2 , £ = —pQ/M, and 
A = arccos(sin0cos£). The sum in Eq. @ splits into 
intra- and interbranch excitations: 

S 2 (E, u) = Y^ A ± 6 ( u - 1 K /v) + B8[uj - (n/n)(q T 2A)], 
Q,q 

(6) 

where Q and q are integer multiples of 2ir, 

A± Sg 2 sin 4 0(1 =F cos£/yT + afn 2 sin 2 £) 2 , 

B = Sg 2 sin 2 0tan 2 0sin 2 £/(l + a 2 K 2 sin 2 £), and 
the sign corresponds to the choice of the initial state. 

For unequal masses £ is incommensurate with tt, and 
one can change summation over Q by integration. The 
first term in Eq. ((6]) leads to strong peaks of S at frequen- 
cies which are integer multiples of I-Kyjlej \i, whereas 
the second term gives lower and wider interbranch lobes 
shown separately as green dashed lines. The widening is 



due to the averaging of A when integrating over Q. 

In the case of equal masses the quantization of the 
center-of-mass motion imposes £ = Q/2 = ■nm. Then, 
the two branches correspond to symmetric and antisym- 
metric two-body states, the interbranch excitations are 
not possible {B = 0) and only the symmetric branch 
is sensitive to the variation of the coupling constant 
[A± = sin 4 0(1 =p (—1)"')]. This completely ignores the 
interbranch peaks and strongly overestimates the intra- 
branch ones, contrary to the numerics. Better agree- 
ment is obtained by assuming the continuum uniform 
distribution of Q as in the nonintegrable cases - all red 
solid curves in Fig. Q] are obtained by integrating over 
Q. The proper distribution of Q and the nature of the 
interbranch peaks is beyond the two-body physics. Yet, 
the binary approximation qualitatively explains the low- 
frequency suppression of S in the integrable case. 

Beyond the case 7 ss 0.5 presented in Fig. [1] we calcu- 
late S(e,ui) for 7 varying from ~ 10~ 2 to ~ 10 2 . Com- 
paring with Fig. [TJ for stronger interactions (smaller 7) 
the interbranch lobes turn into sharper peaks consistent 
with the binary approximation (in this case A is al- 
ways close to 7r/2). For larger 7 the peaks smoothen 
and S(e,oj) is well approximated by S(e, ui) except for 
small w, where S tends to a constant in the nonintegrable 
case. We also calculate S(e,uj) for M = m4 K/ m 6Li an d 
M = TO6Li/w40K and see no qualitative deviations from 

Fig.m 

We have also performed an extensive analysis of the 
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Figure 2: (color online) Dynamic structure factor for the Lieb- 
Liniger model (dotted line) and the binary approximation (red 
solid line). Inset shows the same data on log-log scale. 

bosonic Lieb-Liniger model given by Eq. ([1]) with mi = 1 
and g^ = g. We find that here the two-body approxima- 
tion works perfectly well (with Q = 2-Km) for all coupling 
constants and particle numbers (N < 4) that we have 
considered. Figure [2] is a representative plot of S(e,u) 
for a four-boson system. We observe that at low frequen- 
cies S oc lj 4 (see inset in Fig. [2]). This is, actually, a 
manifestation of the ID fcrmionization as small frequen- 
cies correspond to small relative momenta k = [iuj/2-k 
[see Eq. (j6|)] . The probability to find two atoms at small 
distances drops as n 2 a\, which directly transfers to the 
matrix elements of the perturbation and eventually leads 
to the K 4 -behavior of the coefficient A± in Eq. ([6]). 

Apart from this cj 4 -suppression at small frequencies 
the binary approximation prohibits excitations with u> < 
1/L 2 , which is the minimal level spacing for the two- 
body problem. Assuming that this result holds for large 
N oc L, we can conclude that the original state of the 
system is preserved if the perturbation is slow polynomi- 
ally in the system size. This is consistent with the state- 
ment on the absence of adiabaticity in ID systems in the 
thermodynamic limit (20| . Note, however, the distinction 
between polynomially long timescales for integrable sys- 
tems and exponentially long ones in nonintegrable cases, 
which is interesting from the quantum computing per- 
spective. 

Our particle number is limited mostly by the complex- 
ity of the nonintegrable model. There is still some room 
for increase, but we expect no qualitative change of the 
system behavior, at least at small uj where the results are 
consistent with the random matrix theory. In contrast, 
going to larger N in the integrable cases and develop- 
ing a smarter approach for calculating matrix elements 
seems to be an interesting theoretical project. Caux and 
Calabrese have recently proposed an efficient numerical 
algorithm for calculating correlation functions based on 
the algebraic Bethe ansatz (2lj . 
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